
!  Rheology (Update stresses depending on rheology)
!  Calculate total finite strain and plastic strain  
    
subroutine fl_rheol
include 'precision.inc' 
include 'params.inc'
include 'arrays.inc'

dimension depl(4)
dimension s11p(4),s22p(4),s12p(4),s33p(4),s11v(4),s22v(4),s12v(4),s33v(4)
dimension njTinj(nz)
logical rh_sel
real*8 sarc11,sarc22,tinj, finj, finj2

c1d3 = 1./3.
c1d6 = 1./6.
pi = 3.14159265358979323846			!testing

! initialization of random generators
jran = 13

if( mod(nloop,10).eq.0 .OR. ireset.eq.1 ) then
    rh_sel = .true.
else
    rh_sel = .false.
endif
rh_sel = .true.
irh=irheol(mphase)
if(irh.eq.11) call init_visc
if(iynts.eq.1) call init_temp	!G.Ito 9/06		

!---------------------------------------------------------------------------------------
! Axis Accretion G. Ito 10/06
!---------------------------------------------------------------------------------------
if (ny_inject.gt.0) then
  tinj=tinj+dt
  t_Mout=t_Mout+dt
  finj=1.0
  finj2=1.0
  if (rateM_mech .lt. 1.0.and.ny_inject.eq.3) then
    finj=(tinj-ptect)/ptrans
    finj=dmax1(0.0,finj)
    finj=dmin1(1.0,finj)
    finj2=(tinj-(ptect+ptrans))/ptrans
    finj2=dmax1(0.0,finj2)
    finj2=dmin1(1.0,finj2) 
  endif
  ninjtop=jinj1
  ninjbot=jinj2
  ninj=iinj2-iinj1+1
  nelem_inject=jinj2-jinj1+1

  if (ny_inject.gt.0.and.tinj.ge.period_inj) then
    tinj=0.0
  endif
  if (ny_inject.eq.3) then
  if ((rateM_mech.lt.1.0 .and. t_Mout .ge. 0.1*ptrans) &
    .or. (rateM_mech.ge.1.0.and. t_Mout.ge.10.*ptrans)) then				!G.Ito	
    iacc1=iinj1-1
    iacc2=iinj2+1
    xx1=(cord(1,iacc1,1)+cord(2,iacc1,1)+cord(3,iacc1,1))/3.0
    xx2=(cord(1,iacc2,1)+cord(2,iacc2,1)+cord(3,iacc2,1))/3.0
    dxacc1=((xx2-xx1)-dxaccret1)/2.0
    write(6,'(i4, 4E14.6)') 222, time/sec_year/1.e+3, dxacc1, vbc*sec_year, xx1
    dxaccret1=(xx2-xx1)
    t_Mout=0.0
  endif
  endif
endif
if( mod(nloop,2000) .eq. 0 ) then				!G.Ito	
  dum=(cord(1,iinj11,2)+cord(1,iinj22,2))/2.0  
  write(6,'(i4,i12,3F14.6, 1E13.4, 2F12.4, E12.4, F11.3)') &
    111,nloop, time/sec_year/1.e+3, tinj/sec_year/1.e+3, &
    angemint, boffave, finj, finj2, strain_inert0*strain_inert, dum	
endif

!---------------------------------------------------------------------------------------
! Cheat to initiate a single fault
!---------------------------------------------------------------------------------------
!dtapszero=dtapszero+dt
!if (dtapszero.lt.(0.5*ptect)) then
!  do i=1,int(nx/2)	!temporary
!    write(*,*) aps(i,3), dtapszero
!  do j=1,nz-1		!temporary
!    aps(j,i)=0.0	!temporary
!  end do			!temporary
!  end do			!temporary
!endif
!---------------------------------------------------------------------------------------
!  End Accretion code
!---------------------------------------------------------------------------------------

irh_mark = 0
!!!DIR PREFER_PARALLEL
!$DIR LOOP_PARALLEL
!$DIR LOOP_PRIVATE(j,k,irh,iph,stherm,tmpr,bulkm,rmu,vis,coh,phi,psi,ipls,diss,hardn,sII_plas,sII_visc,quad_area,s0,s0a,s0b)
!$DIR LOOP_PRIVATE(de11,de22,de12,de33,dv,s11p,s22p,s12p,s33p,s11v,s22v,s12v,s33v,depl)
do 3 i = 1,nx-1
   rogh = 0. !Reset for pressure calculation [M. Behn April 4, 2006]
   ! Calculate ninjbot to be the minimum of ninjbot and the depth to the maximum cutoff temp (Tinj)
   jcnt=1
   njTinj(1) = ninjbot
    do j = 1,nz-1
        dcord = cord(1,i,2) - cord(j+1,i,2)
         if(dcord.gt.Tcinj) then
           njTinj(jcnt) = j
           jcnt = jcnt+1
         endif
    end do
    ninjbot = min(ninjbot,njTinj(1))

    do 3 j = 1,nz-1
        ! phasez (j,i) is number of a phase NOT a rheology 
        iph = iphase(i,j,phasez(j,i))
        irh = irheol(iph)

        ! Elastic modules & viscosity & plastic properties
        bulkm = rl(iph) + 2.*rm(iph)/3.
        rmu   = rm(iph)
        coh   = coha(iph)
        phi   = phimean(iph)
        psi   = psia(iph)
        poiss = 0.5*rl(iph)/(rl(iph)+rm(iph))
        young = rm(iph)*2.*(1.+poiss) 

        if (ny_inject.gt.0.and.ny_inject.le.2) then      
	if (i.ge.iinj1.and.i.le.iinj2.and.j.le.ninjbot.and.j.ge.ninjtop) then
           dxinj=0.5*(cord(j,i+1,1)-cord(j,i,1)+cord(j+1,i+1,1)-cord(j+1,i,1))
	   fdum=ratfac*rate_inj_mech*dt/(dxinj*dble(ninj))
           sarc1v = -young/(1.-poiss*poiss)*fdum                        ! Plane Stress 
           sarc2v = poiss*sarc1v
           sarc3v = 0.

           sarc1p = -(1.-poiss)*young/(1.+poiss)/(1.-2.*poiss)*fdum     ! Plane Strain Dikes - M. Behn 04/04/05
           sarc2p = sarc1p*poiss/(1.-poiss)
           sarc3p = poiss * (sarc1p + sarc2p)                                           ! Correct s33 - M. Behn 04/19/05     
        endif
        endif

        ! Thermal stresses (alfa_v = 3.e-5 1/K)
        stherm = 0.
        if (istress_therm.gt.0) stherm = -alfa(iph)*bulkm*(temp(j,i)-temp0(j,i))
          ! Preparation of plastic properties
        if (irh.eq.6 .or. irh.ge.11) call pre_plast(i,j,coh,phi,psi,jran,hardn)
        ! Re-evaluate viscosity
        if (irh.eq.3 .or. irh.eq.12) then 
           if( mod(nloop,ifreq_visc).eq.0 .OR. ireset.eq.1 ) visn(j,i) = Eff_visc(i,j)
        endif	
!------------- For ny_inject=3 ------------------------------------
        if (ny_inject.eq.3.and.tinj.gt.ptect) then
	  if (i.ge.iinj1.and.i.le.iinj2.and.j.le.ninjbot.and.j.ge.ninjtop) then
             tmpr = 0.25*(temp(j,i)+temp(j+1,i)+temp(j,i+1)+temp(j+1,i+1))
	     densT = Eff_dens(i,j)			!G. Ito to keep consistent with rest of code
             dh1 = cord (j,i  ,2) - cord (j+1,i  ,2)
             dh2 = cord (j,i+1,2) - cord (j+1,i+1,2)
             dh  = 0.5 * (dh1+dh2)
             dPT = densT * g * dh
	     !dont know what the next line does here and in rest of code but keep it here for consistency G. Ito
             dP = dPT * ( 1 - beta(iph)*rogh ) / ( 1 + beta(iph)/2*dPT )  
             press = rogh + 0.5*dP
             rogh = rogh + dP
             if (visn(j,i).gt.visinj) then
	       visn(j,i)=visn(j,i)*(visinj/visn(j,i))**finj2  !1/19/07
             endif
!	     if (ireset.eq.1.and.i.ge.iinj11.and.i.le.iinj22) aps(j,i)=0.0d0   !Fix aps during accretion 6/4/07
	  endif
        endif	
!------------- END ny_inject=3 ------------------------------------
	vis = visn(j,i)

        ! Cycle by triangles
        do 4 k = 1,4   

            ! Incremental strains
            de11 = strainr(1,k,j,i)*dt
            de22 = strainr(2,k,j,i)*dt
            de12 = strainr(3,k,j,i)*dt
            de33 = 0.
            dv = dvol(k,j,i)
            s11p(k) = stress0(1,k,j,i) + stherm 
            s22p(k) = stress0(2,k,j,i) + stherm 
            s12p(k) = stress0(3,k,j,i) 
            s33p(k) = stress0(4,k,j,i) + stherm 
            s11v(k) = s11p(k)
            s22v(k) = s22p(k)
            s12v(k) = s12p(k)
            s33v(k) = s33p(k)
            if (irh.eq.1) then
                ! elastic
                call elastic(bulkm,rmu,s11p(k),s22p(k),s33p(k),s12p(k),de11,de22,de12)
                irheol_fl(j,i) = 0  
                stress0(1,k,j,i) = s11p(k)
                stress0(2,k,j,i) = s22p(k)
                stress0(3,k,j,i) = s12p(k)
                stress0(4,k,j,i) = s33p(k)

            elseif (irh.eq.3) then
                ! viscous
                call maxwell(bulkm,rmu,vis,s11v(k),s22v(k),s33v(k),s12v(k),de11,de22,de33,de12,dv)
                irheol_fl(j,i) = -1  
                stress0(1,k,j,i) = s11v(k)
                stress0(2,k,j,i) = s22v(k)
                stress0(3,k,j,i) = s12v(k)
                stress0(4,k,j,i) = s33v(k)

            elseif (irh.eq.6) then
                ! plastic
                call plastic(bulkm,rmu,coh,phi,psi,depl(k),ipls,diss,hardn,&
                     s11p(k),s22p(k),s33p(k),s12p(k),de11,de22,de33,de12)
                irheol_fl(j,i) = 1
                stress0(1,k,j,i) = s11p(k)
                stress0(2,k,j,i) = s22p(k)
                stress0(3,k,j,i) = s12p(k)
                stress0(4,k,j,i) = s33p(k)

            elseif (irh.ge.11) then 
                ! Mixed rheology (Maxwell or plastic)
                if( rh_sel ) then
                    call plastic(bulkm,rmu,coh,phi,psi,depl(k),ipls,diss,hardn,&
                        s11p(k),s22p(k),s33p(k),s12p(k),de11,de22,de33,de12)
                    call maxwell(bulkm,rmu,vis,s11v(k),s22v(k),s33v(k),s12v(k),&
                        de11,de22,de33,de12,dv)
                else ! use previously defined rheology
                    if( irheol_fl(j,i) .eq. 1 ) then
                        call plastic(bulkm,rmu,coh,phi,psi,depl(k),ipls,diss,hardn,&
                            s11p(k),s22p(k),s33p(k),s12p(k),de11,de22,de33,de12)
                        stress0(1,k,j,i) = s11p(k)
                        stress0(2,k,j,i) = s22p(k)
                        stress0(3,k,j,i) = s12p(k)
                        stress0(4,k,j,i) = s33p(k)
                    else  ! irheol_fl(j,i) = -1
                        call maxwell(bulkm,rmu,vis,s11v(k),s22v(k),s33v(k),s12v(k),&
                            de11,de22,de33,de12,dv)
                        stress0(1,k,j,i) = s11v(k)
                        stress0(2,k,j,i) = s22v(k)
                        stress0(3,k,j,i) = s12v(k)
                        stress0(4,k,j,i) = s33v(k)
                    endif
                endif
             endif

4       continue

             if( irh.ge.11 .AND. rh_sel ) then
                ! deside - elasto-plastic or viscous deformation
                sII_plas = (s11p(1)+s11p(2)+s11p(3)+s11p(4)-s22p(1)-s22p(2)-s22p(3)-s22p(4))**2 &
                     + 4*(s12p(1)+s12p(2)+s12p(3)+s12p(4))**2
                sII_visc = (s11v(1)+s11v(2)+s11v(3)+s11v(4)-s22v(1)-s22v(2)-s22v(3)-s22v(4))**2 &
                     + 4*(s12v(1)+s12v(2)+s12v(3)+s12v(4))**2

                if (sII_plas .lt. sII_visc) then
                   do k = 1, 4
                      stress0(1,k,j,i) = s11p(k)
                      stress0(2,k,j,i) = s22p(k)
                      stress0(3,k,j,i) = s12p(k)
                      stress0(4,k,j,i) = s33p(k)

                      if(ny_inject.gt.0.and.ny_inject.le.2.and.j.le.ninjbot.and.j.ge.ninjtop) then

                         if(i.ge.iinj1.and.i.le.iinj2) then
                            if(ntap.eq.0) then
                               dtpr=1.
                            elseif(ntap.eq.1) then 
                               dtpr = (j + 0.5 - ninjtop) * (2./nelem_inject)
                            elseif(ntap.eq.2) then 
                               dtpr = sqrt( (nelem_inject/2.)**2 - ((j-0.5) - ((ninjtop+ninjbot-1.)/2.))**2 ) / (nelem_inject/2.)
                            endif
                            stress0(1,k,j,i) = s11p(k) +sarc1p*dtpr
                            stress0(2,k,j,i) = s22p(k) +sarc2p*dtpr
                            stress0(4,k,j,i) = s33p(k) +sarc3p*dtpr   ! Reset s33 for accretion M. Behn 04/19/05
                         endif

                      elseif(ny_inject.eq.3.and.tinj.gt.ptect.and.j.le.ninjbot.and.j.ge.ninjtop) then

                         if(i.ge.iinj11.and.i.le.iinj22) then
			    ds11=finj*(s11p(k)+press)
			    ds22=finj*(s22p(k)+press)
			    ds33=finj*(s33p(k)+press)			    
			    !magma relieves non-lithostatic tension but cant be more compressive
			    stress0 (1,k,j,i) = dmin1(s11p(k), s11p(k)-ds11)
                            stress0 (2,k,j,i) = dmin1(s22p(k), s22p(k)-ds22)           
                            stress0 (4,k,j,i) = dmin1(s33p(k), s33p(k)-ds33)
			 endif
                          
                      endif
                   end do
                   irheol_fl (j,i) = 1

                else 

                   do k = 1, 4
                      stress0(1,k,j,i) = s11v(k)
                      stress0(2,k,j,i) = s22v(k)
                      stress0(3,k,j,i) = s12v(k)
                      stress0(4,k,j,i) = s33v(k)

                      if(ny_inject.gt.0.and.ny_inject.le.2.and.j.le.ninjbot.and.j.ge.ninjtop) then
                                               
                         if(i.ge.iinj1.and.i.le.iinj2) then
                            if(ntap.eq.0) then
                               dtpr=1.
                            elseif(ntap.eq.1) then 
                               dtpr = (j + 0.5 - ninjtop) * (2./nelem_inject)
                            elseif(ntap.eq.2) then
                               dtpr = sqrt( (nelem_inject/2.)**2 - ((j-0.5) - ((ninjtop+ninjbot-1.)/2.))**2 ) / (nelem_inject/2.)
                            endif
                            stress0(1,k,j,i) = s11v(k) +sarc1v*dtpr
                            stress0(2,k,j,i) = s22v(k) +sarc2v*dtpr
                            stress0(4,k,j,i) = s33v(k) +sarc3v*dtpr ! Reset s33 for accretion M. Behn 04/19/05
                          endif

                      elseif(ny_inject.eq.3.and.tinj.gt.ptect.and.j.le.ninjbot.and.j.ge.ninjtop) then
                         if(i.ge.iinj11.and.i.le.iinj22) then
			    ds11=finj*(s11v(k)+press)
			    ds22=finj*(s22v(k)+press)
			    ds33=finj*(s33v(k)+press)
			    !magma relieves non-lithostatic tension but cant be more compressive			    
			    stress0 (1,k,j,i) = dmin1(s11v(k), s11v(k)-ds11)
                            stress0 (2,k,j,i) = dmin1(s22v(k), s22v(k)-ds22)           
                            stress0 (4,k,j,i) = dmin1(s33v(k), s33v(k)-ds33)
                         endif
                      endif
                   end do
                   irheol_fl (j,i) = -1
                endif
             endif

        ! Averaging of isotropic stresses for pair of elements
        if (mix_stress .eq. 1 ) then
        
            ! For A and B couple:
            ! area(n,it) is INVERSE of "real" DOUBLE area (=1./det)
            quad_area = 1./(area(1,j,i)+area(2,j,i))
            s0a=0.5*(stress0(1,1,j,i)+stress0(2,1,j,i))
            s0b=0.5*(stress0(1,2,j,i)+stress0(2,2,j,i))
            s0=(s0a*area(2,j,i)+s0b*area(1,j,i))*quad_area
            stress0(1,1,j,i) = stress0(1,1,j,i) - s0a + s0
            stress0(2,1,j,i) = stress0(2,1,j,i) - s0a + s0
            stress0(1,2,j,i) = stress0(1,2,j,i) - s0b + s0
            stress0(2,2,j,i) = stress0(2,2,j,i) - s0b + s0

            ! For C and D couple:
            quad_area = 1./(area(3,j,i)+area(4,j,i))
            s0a=0.5*(stress0(1,3,j,i)+stress0(2,3,j,i))
            s0b=0.5*(stress0(1,4,j,i)+stress0(2,4,j,i))
            s0=(s0a*area(4,j,i)+s0b*area(3,j,i))*quad_area
            stress0(1,3,j,i) = stress0(1,3,j,i) - s0a + s0
            stress0(2,3,j,i) = stress0(2,3,j,i) - s0a + s0
            stress0(1,4,j,i) = stress0(1,4,j,i) - s0b + s0
            stress0(2,4,j,i) = stress0(2,4,j,i) - s0b + s0

        endif
        !  ACCUMULATED PLASTIC STRAIN
        ! Average the strain for pair of the triangles
        ! Note that area (n,it) is inverse of double area !!!!!
        if (sII_plas .lt. sII_visc) then			!G.Ito 12/10/06, added if statement
          aps(j,i) = aps(j,i) &
            + 0.5*( depl(1)*area(2,j,i)+depl(2)*area(1,j,i) ) / (area(1,j,i)+area(2,j,i)) &
            + 0.5*( depl(3)*area(4,j,i)+depl(4)*area(3,j,i) ) / (area(3,j,i)+area(4,j,i))
          if( aps(j,i) .lt. 0. ) aps(j,i) = 0.
	endif
        ! LINEAR HEALING OF THE PLASTIC STRAIN
        if (tau_heal .ne. 0.) &
            aps (j,i) = aps (j,i)/(1.+dt/tau_heal) 
        ! TOTAL FINITE STRAIN  
        strain (1,j,i) = strain(1,j,i) + 0.25*dt*(strainr(1,1,j,i)+strainr(1,2,j,i)+strainr(1,3,j,i)+strainr(1,4,j,i))
        strain (2,j,i) = strain(2,j,i) + 0.25*dt*(strainr(2,1,j,i)+strainr(2,2,j,i)+strainr(2,3,j,i)+strainr(2,4,j,i))
        strain (3,j,i) = strain(3,j,i) + 0.25*dt*(strainr(3,1,j,i)+strainr(3,2,j,i)+strainr(3,3,j,i)+strainr(3,4,j,i))
3 continue

return
end
